#include <iostream>
#include <cstring>
#include <map>
#include "TH1D.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "DelphesClasses.h"
#include "ExRootTreeReader.h"
#include "TClonesArray.h"

#include "CmdLine/CmdLine.hh"
using namespace std;


int main(int argc, char** argv)
{
  //initialize input
  CmdLine cmdline(argc, argv);

  //get input and output
  string input, output;
  
  try
    {
      input = cmdline.value<string> ("-i");
      output = cmdline.value<string> ("-o");
    }
  catch(...)
    {
      cout<<"ERROR: please specify input and output via -i input -o output"<<endl;
      return 0;
    }

  //initialize plots
  map<string, TH1D*> hist;

  //book histograms
  hist["nj100"]=new TH1D("nj100", "Number of Jets with p_{T} > 100 GeV", 20, 0, 20);
  hist["nj100_after"]=new TH1D("nj100_after", "Number of Jets with p_{T} > 100 GeV", 20, 0, 20);

  hist["nbj100"]=new TH1D("nbj100", "Number of B Jets with p_{T} > 100 GeV", 10, 0, 10);
  hist["nbj100_after"]=new TH1D("nbj100_after", "Number of B Jets with p_{T} > 100 GeV", 20, 0, 20);


  hist["nj150"]=new TH1D("nj150", "Number of Jets with p_{T} > 150 GeV", 20, 0, 20);

  hist["HT"]=new TH1D("HT", "H_{T}", 10000/200, 0, 10000);
  hist["HT_after"]=new TH1D("HT_after", "H_{T} (after cut)", 10000/200, 0, 10000);
  hist["MJ"]=new TH1D("MJ", "M_{J}", 10000/200, 0, 10000);

  hist["met"]=new TH1D("met", "#slash{E}_{T}", 10000/200, 0, 10000);
  hist["metsig"]=new TH1D("metsig", "#slash{E}_{T}/#sqrt{H_{T}}", 200, 0, 150);

  hist["met_after"]=new TH1D("met_after", "#slash{E}_{T}", 10000/200, 0, 10000);
  hist["metsig_after"]=new TH1D("metsig_after", "#slash{E}_{T}/#sqrt{H_{T}} (after cut)", 200, 0, 150);


  hist["jet1pt"]=new TH1D("jet1pt", "1^{st} jet p_{T}", 10000/200, 0, 10000);
  hist["jet1pt_after"]=new TH1D("jet1pt_after", "1^{st} jet p_{T} (after cut)", 10000/200, 0, 10000);
  hist["jet2pt"]=new TH1D("jet2pt", "2^{nd} jet p_{T}", 10000/200, 0, 10000);
  hist["jet2pt_after"]=new TH1D("jet2pt_after", "2^{nd} jet p_{T} (after cut)", 10000/200, 0, 10000);
  hist["jet3pt"]=new TH1D("jet3pt", "3^{rd} jet p_{T}", 10000/200, 0, 10000);
  hist["jet4pt"]=new TH1D("jet4pt", "4^{th} jet p_{T}", 10000/200, 0, 10000);
  
  hist["dphi"]=new TH1D("dphi", "min #Delta #phi between #slash{E}_{T} and all jet", 100, 0, 3.1415926);

  hist["dphi_after"]=new TH1D("dphi_after", "min #Delta #phi between #slash{E}_{T} and all jet (after cut)", 100, 0, 3.1415926);
  


  //create a TChain
  TChain c("Delphes");
  c.AddFile(input.c_str());
  
  //now initiate an ExRootTreeReader for reading
  ExRootTreeReader* reader=new ExRootTreeReader(&c);
  
  //get some trees
  const TClonesArray *cjets= reader->UseBranch("Jet");
  const TClonesArray *cmet= reader->UseBranch("MissingET");
  const TClonesArray *cevt= reader->UseBranch("Event");

  double weight_scale=1.0/reader->GetEntries();
  
  double min_pt=100;

  for(Long64_t i=0; i< reader->GetEntries(); i++)
    {
      reader->ReadEntry(i);
      
      //get weights
      TIter itevt((TCollection*)cevt);
      LHEFEvent* evt= (LHEFEvent*) itevt.Next();
      double weight=evt->Weight * weight_scale; 

      //initialize variables
      int nj100=0;
      int nj150=0;
      //number of bjets
      int nbj100=0;

      double HT=0;
      double MJ=0;

      //missing Et
      TIter itmet((TCollection*)cmet);
      MissingET* pmet= (MissingET*) itmet.Next();

      double met=pmet->MET;
      TLorentzVector vmet;
      vmet.SetPtEtaPhiM(pmet->MET, 0, pmet->Phi, 0);
      

      double metsig;

      double dphi=10;
      //leading and second jet
      TLorentzVector j1, j2, j3, j4;
      

      //iterate over jets
      TIter itjets((TCollection*)cjets); Jet* jet;
      for(int ijet=0; (jet= (Jet*) itjets.Next()); ijet++)
	if(jet->PT > 20)
	  {
	    TLorentzVector j;
	    j.SetPtEtaPhiM(jet->PT, jet->Eta, jet->Phi, jet->Mass);
	    
	    HT+=j.Mt();
	      
	    if(jet->PT > min_pt )//&& fabs(jet->Eta) < 3.0)
	      {
		nj100++;
		
		if(jet->BTag)
		  nbj100++;

		MJ+=jet->Mass;
	
		double cphi=fabs(j.DeltaPhi(vmet));
		if(cphi < dphi)
		  dphi=cphi;
	
		switch(nj100)
		  {
		  case 1: 
		    hist["jet1pt"]->Fill(jet->PT, weight);
		    j1=j;
		    break;
		      
		  case 2: 
		    hist["jet2pt"]->Fill(jet->PT, weight);
		    j2=j;
		    break;
		    
		  case 3: 
		    hist["jet3pt"]->Fill(jet->PT, weight);
		    j3=j;
		    break;
		    
		  case 4: 
		    hist["jet4pt"]->Fill(jet->PT, weight);
		    j4=j;
		    break;
		  }
		
		if(jet->PT > 150)
		  nj150++;
	      }
	  }


      if(nj100 <2) continue;

      metsig=met/sqrt(HT);


      //fill histograms
      hist["HT"]->Fill(HT, weight);
      hist["MJ"]->Fill(MJ, weight);
      hist["nj100"]->Fill(nj100, weight);
      hist["nbj100"]->Fill(nbj100, weight);
      hist["nj150"]->Fill(nj150, weight);
      hist["met"]->Fill(met, weight);
      hist["metsig"]->Fill(metsig, weight);
      hist["dphi"]->Fill(dphi, weight);
      
      
      //cuts
      if(nj100 >= 2 &&
	 nbj100 >= 2 &&
	 j2.M() >= 160 &&
	 dphi >= 0.32 &&
	 met >= 5700 &&
	 metsig >= 66)
	{
	  hist["HT_after"]->Fill(HT, weight);
	  
	  hist["nj100_after"]->Fill(nj100, weight);
	  hist["nbj100_after"]->Fill(nbj100, weight);
	  
	  hist["met_after"]->Fill(met, weight);
	  hist["metsig_after"]->Fill(metsig, weight);
	  hist["dphi_after"]->Fill(dphi, weight);

	  hist["jet1pt_after"]->Fill(j1.Pt(), weight);	  
	  hist["jet2pt_after"]->Fill(j2.Pt(), weight);	  
	}


    }

  

  TFile fout(output.c_str(), "RECREATE");
  map<string, TH1D*>::iterator it=hist.begin();
  for(;it!=hist.end(); ++it)
    {
      it->second->Write();
    }

  fout.Close();
  return 1;
}
